% SimGrievanceModelLoopFinal.m  - many iterations of the numerical Greivance model for generating Greivance Model Part of Table 1
% Guardado and Pennings (2023) "Shock Persistence and the Study of Armed Conflict: Empirical Biases and Some Remedies" International Interactions
% Contact: Steven Pennings (spennings@worldbank.org; steven.pennings@gmail.com)
%       Replicates Greivance model data for Table 1; output Table1data_Greivance.xlsx, which used in ExcelGraphsTables.xlsx (which actually produces the Table)
%       Calls GrievanceLinearAR1.mod GrievanceLinearSeasonal.mod
%       Uses Matlab R2018a (Optimization Toolbox; Statistics & Machine Learning toolbox); Dynare 5.3 (www.dynare.org)



%% Model set up %%
clear
clc
cd C:\Users\wb464833\Dropbox\SeasonalViolence\TheoryNote\Replication\GrievanceModel
addpath C:\dynare\5.3\matlab

% Parameters (we add the "1" as use this to keep variable later on")
sigma=1;        % CRRA risk aversion
sigma1=sigma;
beta=0.99;      % Quarterly discount rate (in the manuscript we call this delta, because beta is the coefficient in a regression)
beta1=beta;
wbar=1;         % Steady State Wage Rate (not logged!)
phi=3;           % Opportunity cost is -phi  
phi1=phi;
psi=0.45;       % Utility weight for violence (doesn't affect log-linearization)
psi1=psi;
truecoeff1=-phi;

[vbar fval] = fsolve(@(x) psi*x^(-1/phi)-wbar^(1-sigma)*(1-x)^-sigma,[0.5])     %Solve for x=steady state violence. This is vbar=0.0727. Taken from Online Appendix 2.2.
cbar=wbar*(1-vbar);
vbar1=vbar;  
cbar1=cbar; 
wbar1=wbar;    

%% Memorandum items - show simple regression of violence on wages generates true opportunity cost parameter (=-3) for seasonal and transient shocks %%
clearvars -except *1
rho1=0;
dynare GrievanceLinearAR1.mod noclearall
[B,BINT] = regress(v,[ones(1000,1) w]);
B_transientreg1=B(2)   % value approximately equal -3


clearvars -except *1
dynare GrievanceLinearSeasonal.mod noclearall
[B,BINT] = regress(v_wSVshock,[ones(40,1) w_wSVshock]);                 % For seasonality, do IRF rather than simulation because seasonal factors are non-random
B_SVreg1=B(2)  % value approximately equal -3

%% Replicate  Table 1 (Grievance model part; selected commodities) 
rhoV_annual_commodities=[0.863671000000000;0.956483000000000;0.899903000000000;0.868851000000000;0.949999000000000;0.970253000000000;0.922262000000000;0.958187000000000;0.967264000000000;0.969971000000000;];
N=length(rhoV_annual_commodities);
reg_rhosimV_commodities=nan(N,1);

for ii=1:1:N 
rho1=rhoV_annual_commodities(ii).^0.25;
dynare GrievanceLinearAR1.mod noclearall
[B,BINT] = regress(v,[ones(1000,1) w]);
reg_rhosimV_commodities(ii)=B(2);
end

Table1data_Grievance=[rhoV_annual_commodities, reg_rhosimV_commodities];
xlswrite('Table1data_Grievance.xlsx',Table1data_Grievance,'Table1data_Grievance')

